#!/usr/bin/env python3
"""
Spigot Pi Calculator

This script implements various spigot algorithms for calculating π digit by digit.
These algorithms are interesting because they don't use floating-point arithmetic
and can generate arbitrary digits of π without calculating all previous digits.

Algorithms implemented:
1. Unbounded Spigot algorithm (Rabinowitz-Wagon)
2. Gosper's algorithm (continued fraction method)
3. Chudnovsky algorithm (high precision calculation)
4. Digit extraction with visualization

The script also includes performance comparisons and visualizations of the outputs.
"""

import time
import math
import sys
from decimal import Decimal, getcontext


def rabinowitz_wagon_spigot(digits=100):
    """
    Unbounded Spigot algorithm for π calculation.
    
    This algorithm generates digits of π one at a time without
    requiring arbitrary precision arithmetic.
    
    Args:
        digits: Number of digits to compute
    
    Returns:
        String representation of π to the specified number of digits
    """
    # Length of array needed
    length = 10 * digits // 3 + 2
    
    # Initialize array of 2s
    a = [2] * length
    
    # Initialize other variables
    nines = predigit = 0
    result = "3."
    
    for j in range(1, digits + 1):
        # Multiply by 10
        q = 0
        for i in range(length - 1, 0, -1):
            a[i] = a[i] * 10 + q
            q = a[i] // (2 * i - 1)
            a[i] %= (2 * i - 1)
        
        a[1] = a[1] * 10 + q
        q = a[1] // 2
        a[1] %= 2
        
        if q == 9:
            nines += 1
        elif q == 10:
            result += str(predigit + 1) + "9" * nines
            predigit = 0
            nines = 0
        else:
            result += str(predigit) + "9" * nines + str(q)
            predigit = q
            nines = 0
    
    return result


def gosper_continued_fraction(iterations=100):
    """
    Gosper's continued fraction algorithm for π.
    
    This uses the continued fraction representation:
    π = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + ...))))
    
    Args:
        iterations: Number of iterations
    
    Returns:
        Approximation of π
    """
    # Coefficients for continued fraction (pre-computed)
    coeffs = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2]
    
    # Extend with pattern if more iterations needed
    if iterations > len(coeffs):
        pattern = [1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2]
        while len(coeffs) < iterations:
            coeffs.extend(pattern)
    
    # Compute continued fraction from right to left
    value = Decimal(0)
    for i in range(min(iterations, len(coeffs)) - 1, -1, -1):
        value = Decimal(coeffs[i]) + (Decimal(1) / value if value != 0 else 0)
    
    return value


def chudnovsky_algorithm(iterations=10):
    """
    Chudnovsky algorithm for high-precision π calculation.
    
    This is one of the most efficient algorithms for calculating π,
    adding roughly 14 decimal digits per iteration.
    
    Args:
        iterations: Number of terms in the sum
    
    Returns:
        High-precision approximation of π
    """
    getcontext().prec = 100 + 14 * iterations  # Set precision
    
    total = Decimal(0)
    multiplier = Decimal(426880) * Decimal(10005).sqrt()
    
    for k in range(iterations):
        try:
            numerator = Decimal((-1) ** k) * Decimal(math.factorial(6 * k)) * Decimal(13591409 + 545140134 * k)
            denominator = Decimal(math.factorial(3 * k)) * Decimal(math.factorial(k) ** 3) * Decimal(640320 ** (3 * k))
            total += numerator / denominator
        except (OverflowError, ValueError):
            print(f"Overflow at iteration {k}, stopping computation")
            break
    
    return multiplier / total


def pi_digit_extraction(n=100):
    """
    Extract the nth digit of π using the BBP formula.
    
    This implements the Bailey-Borwein-Plouffe formula to extract
    specific hexadecimal digits of π without computing previous digits.
    
    Args:
        n: Position of the hexadecimal digit to extract (0-indexed)
    
    Returns:
        Hexadecimal digit of π at position n
    """
    # Implement hexadecimal digit extraction
    def S(j, n):
        # Calculate the sum for the BBP formula
        total = 0
        for k in range(n+1):
            r = 8*k + j
            total = (total + pow(16, n-k, r) * 16 % r) % 1.0
            
        power = n + 1
        while power <= n + 100:
            r = 8 * power + j
            total += pow(16, n - power, r) / r
            total %= 1.0
            power += 1
            
        return total
    
    # Calculate the sum components
    pi_hex = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) % 1.0
    
    # Convert to hexadecimal
    pi_hex = int(pi_hex * 16)
    return "0123456789ABCDEF"[pi_hex]


def compare_methods():
    """
    Compare different methods for calculating π.
    
    Displays execution time, precision, and digits calculated.
    """
    print("Comparing π calculation methods:")
    print("--------------------------------")
    print(f"Built-in math.pi = {math.pi}")
    
    methods = [
        ("Rabinowitz-Wagon Spigot", lambda: rabinowitz_wagon_spigot(50)),
        ("Gosper Continued Fraction", lambda: gosper_continued_fraction(20)),
        ("Chudnovsky Algorithm", lambda: chudnovsky_algorithm(5))
    ]
    
    results = []
    for name, method in methods:
        print(f"\nCalculating π using {name}...")
        start_time = time.time()
        
        try:
            result = method()
            end_time = time.time()
            
            print(f"Result: {result}")
            print(f"Time: {end_time - start_time:.4f} seconds")
            
            # Calculate error if possible
            try:
                error = abs(float(result) - math.pi)
                print(f"Error: {error:.15e}")
            except (ValueError, TypeError):
                # For string results like the spigot algorithm
                if isinstance(result, str) and result.startswith("3."):
                    pi_str = str(math.pi)
                    matching_digits = 0
                    for i in range(min(len(result), len(pi_str))):
                        if result[i] == pi_str[i]:
                            matching_digits += 1
                        else:
                            break
                    print(f"Matching digits: {matching_digits}")
            
            results.append((name, result, end_time - start_time))
        except Exception as e:
            print(f"Error in {name}: {e}")
    
    return results


def visualize_digits(digits=50):
    """
    Generate a simple visualization of π digits.
    
    Displays the first N digits of π with color-coding
    to highlight patterns.
    
    Args:
        digits: Number of digits to visualize
    """
    # Get digits using Rabinowitz-Wagon algorithm
    pi_string = rabinowitz_wagon_spigot(digits)
    
    print("\nVisualization of π digits:")
    print("--------------------------")
    
    # Display digits with simple ASCII art
    for i, digit in enumerate(pi_string):
        if digit == '.':
            print(".", end="")
        else:
            # Create simple ASCII visualization (higher digits get more '#')
            intensity = min(int(digit) + 1, 10)
            print(digit + "#" * intensity)
    
    # Print frequency of each digit
    frequencies = {}
    for digit in pi_string:
        if digit.isdigit():
            frequencies[digit] = frequencies.get(digit, 0) + 1
    
    print("\nDigit frequencies:")
    for digit in sorted(frequencies.keys()):
        bar = "█" * (frequencies[digit] * 50 // digits)
        print(f"{digit}: {bar} ({frequencies[digit]})")


def extract_specific_digits():
    """
    Extract specific digits of π using BBP formula.
    """
    print("\nExtracting specific hex digits of π:")
    print("-----------------------------------")
    
    # Extract some interesting positions
    positions = [0, 1, 10, 100, 1000]
    for pos in positions:
        try:
            digit = pi_digit_extraction(pos)
            print(f"Hex digit at position {pos}: {digit}")
        except Exception as e:
            print(f"Error extracting digit at position {pos}: {e}")


def main():
    """
    Main function to run the π calculation demonstrations.
    """
    print("Spigot Pi Calculator")
    print("===================")
    
    # Set decimal precision for high-accuracy calculations
    getcontext().prec = 100
    
    print("\nCalculating π with different spigot methods...\n")
    
    try:
        # If we have command line arguments
        if len(sys.argv) > 1:
            if sys.argv[1] == "--spigot":
                digits = int(sys.argv[2]) if len(sys.argv) > 2 else 100
                print(f"Calculating {digits} digits of π with Rabinowitz-Wagon spigot:")
                result = rabinowitz_wagon_spigot(digits)
                print(result)
                
            elif sys.argv[1] == "--chudnovsky":
                iterations = int(sys.argv[2]) if len(sys.argv) > 2 else 5
                print(f"Calculating π with Chudnovsky algorithm ({iterations} iterations):")
                result = chudnovsky_algorithm(iterations)
                print(result)
                
            elif sys.argv[1] == "--visualize":
                digits = int(sys.argv[2]) if len(sys.argv) > 2 else 50
                visualize_digits(digits)
                
            elif sys.argv[1] == "--extract":
                extract_specific_digits()
                
            elif sys.argv[1] == "--compare":
                compare_methods()
                
            else:
                print("Unknown option. Available options:")
                print("  --spigot [digits]    - Calculate π using spigot algorithm")
                print("  --chudnovsky [iters] - Calculate π using Chudnovsky algorithm")
                print("  --visualize [digits] - Visualize π digits")
                print("  --extract            - Extract specific π digits")
                print("  --compare            - Compare different methods")
        else:
            # Default: run all demonstrations
            results = compare_methods()
            visualize_digits(30)
            extract_specific_digits()
            
    except KeyboardInterrupt:
        print("\nCalculation interrupted by user.")
    except Exception as e:
        print(f"\nAn error occurred: {e}")


if __name__ == "__main__":
    main()